student_data <- readRDS("../Data/student_2022.rds")
print(colnames(student_data))

# Your data loading, filtering, and transformation logic as provided:
load(here("data/countrycode.rda"))
readRDS(here("data/student_2022.rds")) -> student_2022

# Load the country names, and join
student_country <- left_join(student_2022, countrycode, by = "country")

student <- student_country %>%
  filter(!is.na(gender), !is.na(math), !is.na(stu_wgt))

# Compute average math scores and gender diff
if (!file.exists("data/math_diff_conf_intervals.rda")) {
  math_diff_df <- student %>%
    group_by(gender, country_name) %>%
    summarise(avg = weighted.mean(math, stu_wgt), .groups = "drop") %>%
    ungroup() %>%
    pivot_wider(country_name, names_from = gender, values_from = avg) %>%
    mutate(diff = female - male, country_name = fct_reorder(country_name, diff))

  # Compute bootstrap samples
  set.seed(2024)
  boot_ests_math <- map_dfr(1:100, ~{
    student %>%
      group_by(country_name, gender) %>%
      sample_n(size = n(), replace = TRUE) %>%
      summarise(avg = weighted.mean(math, stu_wgt), .groups = "drop") %>%
      pivot_wider(country_name, names_from = gender, values_from = avg) %>%
      mutate(diff = female - male, country_name = fct_reorder(country_name, diff)) %>%
      mutate(boot_id = .x)
  })

  # Compute bootstrap confidence intervals
  math_diff_conf_intervals <- boot_ests_math %>%
    group_by(country_name) %>%
    summarise(lower = sort(diff)[5],
              upper = sort(diff)[95],
              .groups = "drop") %>%
    left_join(math_diff_df, by = "country_name") %>%
    mutate(country_name = fct_reorder(country_name, diff)) %>%
    mutate(score_class = factor(case_when(
      lower < 0 & upper <= 0 ~ "boys",
      lower < 0 & upper >= 0 ~ "nodiff",
      lower >= 0 & upper > 0 ~ "girls"),
      levels = c("boys", "nodiff", "girls")))

  save(math_diff_conf_intervals, file = here("data/math_diff_conf_intervals.rda"))
} else {
  load(here("data/math_diff_conf_intervals.rda"))
}

math_plot <- ggplot(math_diff_conf_intervals,
                    aes(x = diff, y = fct_reorder(country_name, diff), col = score_class)) +
  geom_point(size = 3) +                                           # Larger points for better visibility
  geom_errorbar(aes(xmin = lower, xmax = upper), width = 0.5) +    # Error bars with appropriate width
  geom_vline(xintercept = 0, color = "gray", linetype = "dashed") +# Dashed line at x=0 for reference
  scale_colour_manual("", values = c("boys" = "#3288bd",
                                     "nodiff" = "#969696",
                                     "girls" = "#f46d43")) +        # Custom colors for groups
  labs(y = "Country", x = "Difference (Girls - Boys)", title = "Math Gender Gap Analysis") + # Axis titles
  theme_minimal() +                                                # Clean theme
  theme(
    axis.text.y = element_text(size = 8, hjust = 1),               # Smaller y-axis text, aligned left
    axis.text.x = element_text(size = 10),                         # Normal size for x-axis labels
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold"), # Centered title, bold
    plot.margin = margin(20, 20, 20, 150)                          # Add margin for better readability
  ) +
  scale_x_continuous(limits = c(-70, 70),                          # Set x-axis limits
                     breaks = seq(-60, 60, 20),                   # Define breaks
                     labels = abs(seq(-60, 60, 20))) +             # Show absolute values for x-axis
  annotate("text", x = 55, y = 1, label = "Girls", size = 4, color = "red") + # Label for Girls
  annotate("text", x = -55, y = 1, label = "Boys", size = 4, color = "blue")  # Label for Boys

# Print the plot
print(math_plot)

# Your data loading, filtering, and transformation logic as provided:
load(here("data/countrycode.rda"))
readRDS(here("data/student_2022.rds")) -> student_2022

# Load the country names, and join
student_country <- left_join(student_2022, countrycode, by = "country")

# Subset data and drop missing values for reading scores
student <- student_country %>%
  filter(!is.na(gender), !is.na(read), !is.na(stu_wgt))

# Compute average reading scores and gender diff
if (!file.exists("data/read_diff_conf_intervals.rda")) {
  read_diff_df <- student %>%
    group_by(gender, country_name) %>%
    summarise(avg = weighted.mean(read, stu_wgt), .groups = "drop") %>%
    ungroup() %>%
    pivot_wider(country_name, names_from = gender, values_from = avg) %>%
    mutate(diff = female - male, country_name = fct_reorder(country_name, diff))

  # Compute bootstrap samples
  boot_ests_read <- map_dfr(1:100, ~{
    student %>%
      group_by(country_name, gender) %>%
      sample_n(size = n(), replace = TRUE) %>%
      summarise(avg = weighted.mean(read, stu_wgt), .groups = "drop") %>%
      ungroup() %>%
      pivot_wider(country_name, names_from = gender, values_from = avg) %>%
      mutate(diff = female - male, country_name = fct_reorder(country_name, diff)) %>%
      mutate(boot_id = .x)
  })

  # Compute bootstrap confidence intervals
  read_diff_conf_intervals <- boot_ests_read %>%
    group_by(country_name) %>%
    summarise(lower = sort(diff)[5],
              upper = sort(diff)[95],
              .groups = "drop") %>%
    left_join(read_diff_df, by = "country_name") %>%
    mutate(country_name = fct_reorder(country_name, diff)) %>%
    mutate(score_class = factor(case_when(
      lower < 0 & upper <= 0 ~ "boys",
      lower < 0 & upper >= 0 ~ "nodiff",
      lower >= 0 & upper > 0 ~ "girls"),
      levels = c("boys", "nodiff", "girls")))

  save(read_diff_conf_intervals, file = here("data/read_diff_conf_intervals.rda"))
} else {
  load(here("data/read_diff_conf_intervals.rda"))
}

# Plot reading scores with ggplot2
read_plot <- ggplot(read_diff_conf_intervals,
                    aes(diff, country_name, col = score_class)) +
  scale_colour_manual("", values = c("boys" = "#3288bd",
                                     "nodiff" = "#969696",
                                     "girls" = "#f46d43")) +
  geom_point() +
  geom_errorbar(aes(xmin = lower, xmax = upper), width = 0) +
  geom_vline(xintercept = 0, color = "#969696") +
  labs(y = "", x = "", title = "Reading") +
  theme(legend.position = "none") +
  annotate("text", x = 50, y = 1, label = "Girls") +
  annotate("text", x = -50, y = 1, label = "Boys") +
  scale_x_continuous(limits = c(-70, 70),
                     breaks = seq(-60, 60, 20),
                     labels = abs(seq(-60, 60, 20)))

# Convert the ggplot object to an interactive Plotly plot
read_plotly <- ggplotly(read_plot)

# Display the interactive Plotly plot
read_plotly

# Your data loading, filtering, and transformation logic as provided:
load(here("data/countrycode.rda"))
readRDS(here("data/student_2022.rds")) -> student_2022

# Load the country names, and join
student_country <- left_join(student_2022, countrycode, by = "country")

# Subset data and drop missing values for science scores
student <- student_country %>%
  filter(!is.na(gender), !is.na(science), !is.na(stu_wgt))

# Compute average science scores and gender diff
if (!file.exists("data/sci_diff_conf_intervals.rda")) {
  sci_diff_df <- student %>%
    group_by(gender, country_name) %>%
    summarise(avg = weighted.mean(science, stu_wgt), .groups = "drop") %>%
    ungroup() %>%
    pivot_wider(country_name, names_from = gender, values_from = avg) %>%
    mutate(diff = female - male, country_name = fct_reorder(country_name, diff))

  # Compute bootstrap samples
  boot_ests_sci <- map_dfr(1:100, ~{
    student %>%
      group_by(country_name, gender) %>%
      sample_n(size = n(), replace = TRUE) %>%
      summarise(avg = weighted.mean(science, stu_wgt), .groups = "drop") %>%
      ungroup() %>%
      pivot_wider(country_name, names_from = gender, values_from = avg) %>%
      mutate(diff = female - male, country_name = fct_reorder(country_name, diff)) %>%
      mutate(boot_id = .x)
  })

  # Compute bootstrap confidence intervals
  sci_diff_conf_intervals <- boot_ests_sci %>%
    group_by(country_name) %>%
    summarise(lower = sort(diff)[5], upper = sort(diff)[95], .groups = "drop") %>%
    left_join(sci_diff_df, by = "country_name") %>%
    mutate(country_name = fct_reorder(country_name, diff)) %>%
    mutate(score_class = factor(case_when(
      lower < 0 & upper <= 0 ~ "boys",
      lower < 0 & upper >= 0 ~ "nodiff",
      lower >= 0 & upper > 0 ~ "girls"),
      levels = c("boys", "nodiff", "girls")))

  save(sci_diff_conf_intervals, file = here("data/sci_diff_conf_intervals.rda"))
} else {
  load(here("data/sci_diff_conf_intervals.rda"))
}

# Plot science scores with ggplot2
sci_plot <- ggplot(sci_diff_conf_intervals,
                    aes(diff, country_name, col = score_class)) +
  scale_colour_manual("", values = c("boys" = "#3288bd",
                                     "nodiff" = "#969696",
                                     "girls" = "#f46d43")) +
  geom_point() +
  geom_errorbar(aes(xmin = lower, xmax = upper), width = 0) +
  geom_vline(xintercept = 0, color = "#969696") +
  labs(y = "", x = "", title = "Science") +
  theme(legend.position = "none") +
  annotate("text", x = 50, y = 1, label = "Girls") +
  annotate("text", x = -50, y = 1, label = "Boys") +
  scale_x_continuous(limits = c(-70, 70),
                     breaks = seq(-60, 60, 20),
                     labels = abs(seq(-60, 60, 20)))

# Convert the ggplot object to an interactive Plotly plot
sci_plotly <- ggplotly(sci_plot)

# Display the interactive Plotly plot
sci_plotly

theme_map <- function(...) {
  theme_minimal() +
  theme(
    text = element_text(),
    axis.line = element_blank(),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    panel.border = element_blank(),
    panel.grid = element_blank()
  )
}

region2country = function(region_name){
  country_name = case_when(
    region_name == "Brunei Darussalam" ~ "Brunei",
    region_name == "United Kingdom" ~ "UK",
    region_name %in% c("Macau SAR China", "B-S-J-Z (China)",
                        "Hong Kong SAR China") ~ "China",
    region_name == "Korea" ~ "South Korea",
    region_name == "North Macedonia" ~ "Macedonia",
    region_name == "Baku (Azerbaijan)" ~ "Baku",
    region_name %in% c("Moscow Region (RUS)", "Tatarstan (RUS)",
                        "Russian Federation") ~ "Russia",
    region_name == "Slovak Republic" ~ "Slovakia",
    region_name == "Chinese Taipei" ~ "Taiwan",
    region_name == "United States" ~ "USA",
    TRUE ~ as.character(region_name))
}

math_map_data <- math_diff_conf_intervals  %>%
  dplyr::mutate(country_name = region2country(region_name = country_name)) 

world_map <- map_data("world") %>%
  filter(region != "Antarctica") %>%
  fortify() %>%
  rename(country_name = region)

math_world_data <- full_join(
  x = math_map_data,
  y = world_map,
  by = "country_name") %>% 
  rename(Country = country_name,
         math = diff) %>%
  mutate(math = round(math, digits = 2))

read_map_data <- read_diff_conf_intervals %>%
  dplyr::mutate(country = region2country(region_name = country_name))

world_map <- map_data("world") %>%
  filter(region != "Antarctica") %>%
  fortify() %>%
  rename(country_name = region)

read_world_data <- full_join(
  x = read_map_data,
  y = world_map,
  by = "country_name") %>% 
  rename(Country = country_name,
         Reading = diff) %>%
  mutate(Reading = round(Reading, digits = 2))

sci_map_data <- sci_diff_conf_intervals %>%
  dplyr::mutate(country = region2country(region_name = country_name))

world_map <- map_data("world") %>%
  filter(region != "Antarctica") %>%
  fortify() %>%
  rename(country_name = region)

sci_world_data <- full_join(
  x = sci_map_data,
  y = world_map,
  by = "country_name") %>% 
  rename(Country = country_name,
         Science = diff)  %>%
  mutate(Science = round(Science, digits = 2))

math_dat <- math_world_data %>%
  dplyr::select(Country, math, lat, long, group)

read_dat <- read_world_data %>%
  dplyr::select(Country, Reading, lat, long, group)

sci_dat <- sci_world_data %>%
  dplyr::select(Country, Science, lat, long, group)

math_read_dat <- left_join(math_dat,
                           read_dat,
                           by = c("Country","lat", "long", "group"))

math_read_sci_dat <- left_join(math_read_dat,
                           sci_dat,
                           by = c("Country","lat", "long", "group"))

math_read_sci_dat_wider <- math_read_sci_dat %>%
    pivot_longer(cols = c(2,6,7), names_to = "subjects")

mrs_maps <- ggplot(math_read_sci_dat_wider,
       aes(x = long,
           y = lat,
           group = group)) +
  geom_polygon(aes(fill= value,
                   label = Country)) +
  facet_wrap(~subjects, scales = "free", nrow = 3) +
  theme_map() +
  labs(title = "World Map displaying Gender Gap Scores in Math, Reading and Science")  +
  scale_fill_distiller(palette = "Spectral")

ggplotly(mrs_maps, width=800, height=600)

p1 <- ggplot(data = student_country,
             aes(x = math, y = read)) +
  geom_hex() +
  labs(x = "Math Scores",
       y = "Reading Scores") +
  theme(legend.position="none")

p2 <- ggplot(data = student_country,
             aes(x = math, y = science)) +
  geom_hex() +
  labs(x = "Math Scores",
       y = "Science Scores") +
  theme(legend.position = "none")

p3 <- ggplot(data = student_country,
             aes(x = read, y = science)) +
  geom_hex() +
  labs(x = "Reading Scores",
       y = "Science Scores") +
  theme(legend.position="none")

p1+p2+p3

# Load student data, and filter to country, cache a copy of the data
# to save downloading every time paper is knitted

# if (!file.exists("data/student_all.rda")) {
#   student_all <- load_student("all")
#   save(student_all, file="data/student_all.rda")
# } else {
#   load("data/student_all.rda")
# }
student_all <- load(here("data/student.rda"))
# Give countries their name, subset to four, and select only variables needed
student_country <- left_join(student,
                             countrycode,
                             by = "country") %>%
  dplyr::filter(country_name %in%
                  c("Australia",
                    "Germany",
                    "Peru",
                    "Qatar",
                    "Belgium",
                    "Brazil",
                    "Denmark",
                    "Greece",
                    "Thailand",
                    "Singapore",
                    "Canada",
                    "Portugal")) %>%
  dplyr::select(year, country_name, math, read, science, stu_wgt) %>%
  na.omit() %>%
  pivot_longer(c(math, read, science), names_to = "subject", values_to = "score")

# Compute the bootstrap confidence intervals, and cache result
if (!file.exists("data/all_bs_cf.rda")) {
  all_bootstrap <- map_dfr(1:100, ~{
    student_country %>%
    group_by(country_name, #year,
             subject) %>%
    #sample_n(size = n(), replace = TRUE) %>%
    mutate(year = sample(year, replace=FALSE)) %>%
    group_by(country_name, year,
             subject) %>%
    dplyr::summarise(
      avg = weighted.mean(score, w = stu_wgt, na.rm = TRUE), .groups = "drop") %>%
    #ungroup() %>%
    mutate(boot_id = .x)
  })

  all_bootstrap_ci <- all_bootstrap %>%
    group_by(country_name, year,
             subject) %>%
    summarise(
      lower = min(avg), # sort(avg)[5],
      upper = max(avg), #sort(avg)[95],
      .groups = "drop")

  # compute original estimate of average and join
  all_avg <- student_country %>%
    group_by(country_name, year, subject) %>%
    summarise(
      avg = weighted.mean(score,
                          w = stu_wgt, na.rm = TRUE),
      .groups = "drop")

  all_bs_cf <- left_join(all_avg,
                      all_bootstrap_ci,
                      by = c("country_name",
                             "year",
                             "subject"))

  save(all_bs_cf, file= here("data/all_bs_cf.rda"))

} else {
  load(here("data/all_bs_cf.rda"))
}

all_bs_cf <- all_bs_cf %>%
  mutate(year = as.numeric(as.character(year)),
         country_name = factor(country_name))
                 # levels = c("Singapore",
                 #          "Australia",
                 #          "New Zealand",
                 #          "Germany",
                 #          "Qatar",
                 #          "Indonesia"))


country_names_highlight <- c("Australia", 
                             "Germany", 
                             "Peru", 
                             "Qatar", 
                             "Belgium", 
                             "Brazil", 
                             "Denmark", 
                             "Greece",
                             "Thailand", 
                             "Singapore", 
                             "Canada", 
                             "Portugal")

math_all_bs_cf_plot <- all_bs_cf %>% 
  dplyr::filter(subject == "math") %>% 
  ggplot(aes(x = year, 
             y = avg)) +
  geom_point(alpha = 0.45) +
  geom_line(aes(group = country_name)) +
  gghighlight::gghighlight(country_name %in% country_names_highlight) +
  labs(
    title = "Maths",
     x = "",
     y = "Score") 


read_all_bs_cf_plot <- all_bs_cf %>% 
  dplyr::filter(subject == "read") %>% 
  ggplot(aes(x = year, 
      y = avg)) +
  geom_point(alpha = 0.45) +
  geom_line(aes(group = country_name)) +
  gghighlight::gghighlight(country_name %in% country_names_highlight) +
  labs(
    title = "Reading",
     x = "",
     y = "Score") 

sci_all_bs_cf_plot <- all_bs_cf %>% 
  dplyr::filter(subject == "science") %>% 
  ggplot(aes(x = year, 
      y = avg)) +
  geom_point(alpha = 0.45) +
  geom_line(aes(group = country_name)) +
  gghighlight::gghighlight(country_name %in% country_names_highlight) +
  labs(
    title = "Science",
     x = "",
     y = "Score")


math_all_bs_cf_plot + read_all_bs_cf_plot + sci_all_bs_cf_plot